Coexistence and competition of nematic and gapped states in bilayer graphene 
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In bilayer graphene, the phase diagram in the plane of a strain-induced bare nematic term, No, 
and a top-bottom gates voltage imbalance, Uo, is obtained by solving the gap equation in the 
random-phase approximation. At nonzero No and Uo, the phase diagram consists of two hybrid 
spin-valley symmetry-broken phases with both nontrivial nematic and mass-type order parameters. 
The corresponding phases are separated by a critical line of first- and second-order phase transitions 
at small and large values of No, respectively. The existence of a critical end point, where the line 
of first-order phase transitions terminates, is predicted. For No = 0, a pure gapped state with a 
broken spin- valley symmetry is the ground state of the system. As No increases, the nematic order 
parameter increases, and the gap weakens in the hybrid state. For Uo = 0, a quantum second-order 
phase transition from the hybrid state into a pure gapless nematic state occurs when the strain 
reaches a critical value. A nonzero Uo suppresses the critical value of the strain. The relevance of 
these results to recent experiments is briefly discussed. 



I. INTRODUCTION 

At present, a significant amount of attention is be- 
ing paid to the ground state of bilayer graphene at the 
neutrality point. Various gapped states with broken 
spin-valley symmetry, such as quantum anomalous Hall 
(QAH) ji quantum spin Hall (QSH)^ quantum valley Hall 
(QVH)^ layer antiferromagnet (LAF)^ as well as a gap- 
less nematic stated were suggested as candidates for 
the ground state. Experiments^— showed that bilayer 
graphene at the neutrality point is gapped in the absence 
of external fields. However, the experiment performed by 
the Manchester groupii found a gapless state, and there 
were strong indications that it was a nematic state. 

A variety of suggested ground states is related to the 
way in which the approximate SU(A) spin- valley symme- 
try of the low energy effective model of bilayer graphene 
is broken in its ground state. The application of exter- 
nal electric and magnetic fields adds more complexity 
to the problem of the ground state of bilayer graphene. 
However, it opens the possibility to use them as use- 
ful probes of the ground state of the system. Theoreti- 
cal papers 1 ^— predicted and experiments^ ' 10 ' 15 showed 
the presence of a phase transition between the QSH and 
the QVH states as an electric field increases. Since the 
QVH (layer polarized) state is the ground state of the 
system for a sufficiently large electric field, this rules out 
the QVH state as a candidate for the ground state in the 
absence of external fields. 

Simple continuity arguments suggest that the LAF and 
QSH states are the most likely candidates for the gapped 
ground state of bilayer graphene at the neutrality point 
in the absence of external fields. Indeed, for a sufficiently 
large magnetic field, the QSH (spin polarized) state is the 
ground state of bilayer graphene. On the other hand, 



the experiments^— do not show any phase transitions 
in bilayer graphene at the neutrality point as a magnetic 
field is switched off. This excludes the QAH state because 
it cannot be smoothly connected to the spin polarized 
state at large magnetic fields. The LAF state, on the 
other hand, can adiabatically evolve into the QSH onei& 
and, therefore, cannot be excluded. 

In Refs. [H and [I?], it was shown that, unlike single- 
layer graphene with its linear dispersion relation, the 
quadratic dispersion relation in bilayer graphene leads 
to an instability of its normal phase for arbitrary weak 
Coulomb or short-range repulsive interactions. Rcnor- 
malization group studies in the normal phase of bilayer 
graphene in Refs. H and @ revealed a variety of instabili- 
ties in this phase, which can potentially lead to different 
(competing) ground states with condensates. It was also 
found that the instability with respect to the generation 
of a nematic order parameter is strongest, which seems 
to suggest that the nematic state is the ground state of 
bilayer graphene. However, since the condensates essen- 
tially modify the gap equations, in order to determine the 
phase diagram of the system, it is necessary to analyze 
the gap equations for states with different condensates 
and then to compare their energy densities. This is the 
main goal of the present paper. 

Solving the gap equation in the random-phase approx- 
imation, we find the phase diagram in the plane of a 
strain-induced bare nematic term Ao and a top-bottom 
gates voltage imbalance Uq. We show that the ground 
state of the system at No ^ and C/ ^ is a hybrid 
state with nonzero nematic and mass-type order param- 
eters. The critical line separates the phase with one of 
the hybrid QAH, QSH, or LAF states (which are degen- 
erate in energy in the model at hand) from the hybrid 
QVH state. The phase transition along a large part of 
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this critical line is of second-order except for the region 
of small Nq where a sufficiently large value of U drives 
the first-order phase transition between the two different 
states with broken spin- valley symmetry The predicted 
existence of a critical end point in the phase diagram may 
be relevant to current experiments in bilayer graphene. 

The paper is organized as follows. The analysis is per- 
formed in the framework of the two-band model of bilayer 
graphene with the Coulomb interaction between quasi- 
particles described in Sec. HU In Sec. IIIII we derive the 
gap equations in the random-phase approximation for the 
gapped and nematic order parameters. The main results 
of the paper, including the numerical solutions to the 
gap equations, are presented in Sec. flV] Sec. [V] contains 
a brief summary of the results, as well as our discussions 
and conclusions. In the Appendix the derivation of the 
expression for the energy density is given. 



II. MODEL 

The free part of the effective low-energy Hamiltonian 
of bilayer graphene is as follows J£ 

where n = k\ + ik 2 is the canonical momentum operator, 
m = ji/2vp, where the Fermi velocity is vf — c/300, 
and 71 « 0.4 eV. The two-component spinor field 
carries the valley (£ = ± for K and K' valleys, re- 
spectively) and spin (s = ±) indices. We will use the 
standard convention^ \I/+ S = (V>+Ai , ip+B 2 )s, whereas, 
^-s = (^-Bjj ip-A^s- Here, A\ and B 2 correspond to 
those sublattices in layers 1 and 2, respectively, which, 
according to Bernal (A 2 — B\) stacking, are relevant for 
the low energy dynamics. Let us emphasize that the sub- 
lattice and layer degrees of freedom are not independent 
in this low-energy model: The sublattices A\ and B 2 
correspond to layers 1 and 2, respectively. The effective 
Hamiltonian ((T|) is valid up to energies A = 7i/4«0.1 eV 
and we ignore small trigonal warping effects. 18 

The Coulomb interactions in bilayer graphene are 
given by 3 -^ 

Hc = \ jd 2 xd 2 x'[v{x - x 1 ) [ Pl (x) Pl (x') + p 2 (x)p 2 (x')} 
+2V 12 (x-x') Pl {x)p 2 (x')}, (2) 

where the interaction potentials V{x) and V\ 2 {x) de- 
scribe the intralayer and interlayer interactions, respec- 
tively. Their Fourier transforms are V(k) = 2ne 2 / 'nk 
and Vi 2 (k) = 2ire 2 e~ kd / nk, where d rj 0.35 nm is the 
distance between the layers, and k is a dielectric con- 
stant (in the numerical analysis below, we set k = 4). 
The two-dimensional charge-density operators p± (x) and 



p 2 (x) are 

Pit?) = E *U*Wt)*U*l (I -1,2), (3) 

where P^g) = (1 + £t 3 )/2 and P 2 (g) = (1 - £r 3 )/2 arc 
projectors on states in the layers 1 and 2, respectively. 
Here the matrix £ acts as £^f a = £,^£s, and r 3 is the 
diagonal Pauli matrix acting on the two components of 
the fields Vl'+s and *&- s - Note that the presence of £ in 
-Pi(£) and P 2 (g) is related to the opposite order of the 
A\ and B 2 components in and ty- s . 

Whereas, the gaps in the QAH and QSH states in 
bilayer graphene arc analogous to a Haldane mass^ 
the gaps in the QVH and LAF states are analo- 
gous to a Dirac mass 4^-2. In bilayer graphene, the 
latter is realized as a voltage imbalance between the 
layers ! 12 ' 18 The corresponding order parameters (conden- 
sates) are as follows^^: spin singlets ( , F''T 3 4 , ) (QAH) 
and ("F^t 3 *) (QVH), and spin triplets (^T 3 a 3 ^) 
(QSH) and (*+£r 3 cr 3 *) (LAF), where a 3 is a spin Pauli 
matrix (the indices £ and s in the field are omitted 
here). 

Without a magnetic field, the QSH and QVH states 
are invariant under time reversal,^ whereas, the QAH 
and LAF states break this symmetry^ (the QAH state is 
associated with the Chern number— and the QSH state is 
a two-dimensional topological insulator associated with 
the Z 2 topological invariant). On the other hand, al- 
though the QAH state is invariant under the SU(A), the 
QSH one breaks the spin SU(2) subgroup of SU(A) down 
to U(1)M2. Both the QVH and the LAF states break a 
Z 2 subgroup of the spin- valley Sf/(4) describing the val- 
ley transformation £ — > — £ ("f+s — > ^-s),— m addition 
to that, the LAF state breaks the spin SU(2) down to 
U(l). 

A usual nematic order parameter breaks the SO (2) 
rotational group down to a discrete Z 2 subgroup. As 
we discuss below, however, nematic ordering in bilayer 
graphene may become hybridized with the ordering of 
spin-valley symmetry-broken states. As a result, there 
can exist four different types of nematic ordering. The 
corresponding expressions for the nematic order param- 
eters (condensates) can be formally obtained from those 
in the spin- valley symmetry-broken gapped phases by re- 
placing the matrix t 3 with r , i.e., they are given by 
the following anisotropic order parameters: (^t 1 ^), 
(ttrVf), (^It 1 *), and (^gr 1 ^). These gener- 
alized nematic order parameters break the continuous 
spin-valley symmetry in the same fashion as the corre- 
sponding order parameters in the QAH, QSH, QVH and 
LAF phases, respectively. Also, although two of them, 
(\f r 'l'T 1 er 3 \f r ) and (^^r 1 ^/), are invariant under time re- 
versal, the other two, (^t 1 ^) and (^t^W 3 *) break 
this symmetry. From a physics viewpoint, different types 
of nematic ordering arc possible when quasiparticlcs with 
different spins and valleys contribute unequally to the 
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nematic condensates. In states with broken spin- valley 
symmetries, this should be, of course, generally expected. 



III. QUASIPARTICLE PROPAGATOR AND 
GAP EQUATION 

Our goal is to solve the Schwingcr-Dyson (gap) equa- 
tion for the quasiparticlc propagator in the spin-valley 
symmetry-broken gapped states as well as gapless ne- 



matic states using a unified framework. In general, the 
full quasiparticlc propagator in the coordinate space at 
fixed valley and spin is defined by 

G es (x - x') = -i(Q\T* ts {x)*\ s {x')\Q). (4) 

In the random-phase approximation, the corresponding 
gap equation is readily obtained from the Baym-Kadanoff 
functional in the two-loop approximation, see Eq. (9) in 
the second paper in Ref.ha 



J 



G-^x-x') 



S7}(x - x') - iG is {x - x')V eS (x - x') 



i [Pxi^G^x - x')P 2 {0 + P 2 (Z)Gt s (x - x')Pi(0] V^(x - x') 

l -s*{x~x') [p, (o-p 2 (o\ y, te{[Pi(z')-m')]GM*-^}ViL(p), 



(5) 



r 



where x = (t, x) and the trace is taken over the spinor 
components of the quasiparticle propagator. The dy- 
namically screened interactions V e s and VJl are defined 
by Eqs. (A5) and (A6) in the second paper of Ref. [H, 
Vil(O) = — 2ire 2 d/K is the Fourier transform of the bare 
interaction V^ re (x) taken at u = k = 0. The explicit 
form of the inverse free propagator in momentum space 
reads 



2m 



£U g-e- 

2m 

Af u 




(6) 



where ifk is a polar angle of the quasiparticle momen- 
tum with respect to the orientation of the nematic order, 
Uq = eEj_d/2 is the top-bottom gates voltage imbalance, 
E± is the electric field perpendicular to the layers, and 
A/"q is the bare nematic order parameter due to strain^ 
and rotational mismatch 2 ^ between the layers of bilayer 
graphene. Its value is related, for example, to the angle 
of rotational mismatch 9 as follows^ 



fhv F \ 2 e 2 



V 2a 



7i 



(7) 



where a = 0.142 nm is the intralayer distance between 
neighboring carbon atoms. 

The Fock contribution is given by the second and third 
terms on the right-hand side of Eq. ([5]), whereas, the 
fourth term describes the Hartree contribution. [Note 
that, in accordance with Gauss's law, we omitted the 
Hartree contribution connected with the total electron 
charge of the system proportional to s tr[G^ s (0)], 
which is neutralized by the positive charges of ions 
in the system.] Taking into account that (pi) = 
—i s tr [P(C) G^ s (0)] defines the electron charge den- 
sity in the Ith. layer, we see that the Hartree term propor- 
tional to the Vil(0) interaction describes the layer charge- 
density imbalance. One of the Fock terms in Eq. ([5]) is 
proportional to Vtl- Since Vtl ~ d, it is suppressed com- 
pared to the one with V e g. Because of the polarization 
effects, Vn, is also suppressed compared to p^ are ,2 There- 
fore, it is justifiable to neglect the Fock term with the Vix 
interaction in the analysis that follows. 

Taking into account the gapped and nematic order pa- 
rameters in the inverse quasiparticlc propagator at fixed 
valley and spin, we use the following ansatz: 



G£(w,k) 



k 

2m 



Z-^w.kJw-CA^w.k) |^4( w ,k) e - 2 ^ k) 
'-A(u, k)e 2 ^- + w 6a (u, k) Z- 1 (u, k)w + £A ?S (w, k) 



(8) 



r 



where k and tpk are polar coordinates in the (fci,^) 
plane. Note that although, in general, w^ s is a complex 
function, one can show that, without the loss of general- 
ity, its phase can be set to zero. The functions Z(uj,k) 



and -A(w,k) define quasiparticle residue at the pole and 
rcnormalization of the kinetic term, respectively. 

In Eq. dl]), A^ s = U S +(,A S and W£ s = J\f s +£M S are the 
spin- valley symmetry-breaking and nematic parameters, 
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which arc connected to the order parameters discussed 
above through the following relationship: 

(*+0*) = -*E/ ?j0tr[OGUu,U>)- (9) 

Here the matrix O is t 3 , r 3 cr 3 , £r 3 , and £t 3 ct 3 for the 
QAH, QSH, QVH, and LAF order parameters, respec- 
tively, and O — t 1 , tV 3 , ^t 1 , and ^rV 3 for nematic 
order parameters. [In Eq. ([9]), the trace runs over the val- 
ley degree of freedom.] For the QSH state at Uq = and 
Afo = 0, only the spin-triplet Haldanc mass (A+ = — A_) 
is nonzero in A^ s . On the other hand, for the QAH 
state, only the spin-singlet Haldane mass (A+ = A_) is 
nonzero. Similarly, only the spin-antisymmetric part of 
U s , U + = — £/_, describes the LAF state, and the sym- 
metric in spin voltage U+ = U- occurs in the QVH state. 
In fact, at Uq = these definitions remain unchanged 



even in the corresponding hybrid states, in which an ad- 
ditional spin-singlet nematic parameter A/+ = M- is in- 
duced by the bare nematic term A/o- 

When Uq ^ 0, however, the most general hybrid QSH, 
QAH and LAF phases get further modified. In the hybrid 
QSH state, for example, there will appear a nonzero con- 
tribution of a spin-symmetric parameter [/+ = U- and 
a spin-antisymmetric parameter AA+ = —A4-. In the 
hybrid QAH state, the admixture of two spin-singlet pa- 
rameters U+ = U- and A4 + = A4- will appear. Finally, 
in the hybrid LAF state, the new nonzero contributions 
will be a spin-symmetric parameter [/+ = L/_ and a spin- 
antisymmetric parameter A/+ = — AT- . 

Although A^s and w^ s can, in general, be energy and 
momentum dependent, in our analysis, we use the con- 
stant approximation for them. Also, we set Z(uj,\&) = 
A(u},k) = 1 in G7*(u),K) in Eq.©. Then, the propaga- 
tor reads 



The poles of this propagator determine the dispersion 
relations of quasiparticlcs, 



A- 2 



k 2 



w is (1 + cos2(p k ) + A| s . 

(11) 

As we see, a nonzero dynamical parameter A^ s results in 
a fully gapped dispersion relation for the quasiparticles 
with fixed £ and s. When A^ s = 0, in contrast, the 
quasiparticles with k = y / 2m|w^ s | have a vanishing gap 
for ifk = ±7r/2 and w^ s > 0, and for ip^ = 0, ir when 
w^ s < 0. 

Multiplying the Schwingcr-Dyson equation ([5]) by T3 
and Ti Pauli matrices, respectively, taking the trace over 
the spinor components, and making the Wick rotation, 
we obtain the following set of gap equations for A^ s and 
u> 5s : 



A $a = U - 

W£ a =Af + 
where 



dwd 2 k A is 1 - 

{2ttY D^ s {uj, k) 2 

(12) 



(2^) 3 D f .(w,fc) 

dojd 2 k 2A ?S 



) 3 D is (co,k) 



(14) 



is the charge density imbalance between the layers and 



^ / , s 9 a 9 9 fc 4 k 2 Wf s cos 2ipk 
D ( Mk)=^ + Al+wl + — + ™ . (15) 



In the effective potential V b s(uj, k) in Eqs. (|12p and (|13|). 
we take into account the dynamical polarization of the 
Coulomb interaction^ 



V eff (w, k) = 



2ve 2 



4e 2 m In 4 



(16) 



^l + [Amu ln4/(7rfc 2 )] 2 



At fixed values of t/o and Ao, the gap equations (fT2|) and 
(|13p generally admit many solutions. In order to deter- 
mine which of them corresponds to the ground state, we 
need to find the solution with the lowest energy density. 
The expression for the energy density is obtained in the 
Appendix. 

It should be mentioned that, in the model at hand, 
which includes only the dominant Coulomb interaction 
and has no external magnetic field, the three types of so- 
lutions with spontaneous spin- valley symmetry breaking, 
QAH, LAF, and QSH, appear to be degenerate in energy. 
Moreover, the order parameters in these three phases are 
related by some transformations. This means, in partic- 
ular, that having the explicit form for one of them allows 
for easily reconstructing the other two. For example, if 
the spin-singlet hybrid QAH solution is found to take the 
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following form: 



A+ = A_ = fi(U ,Afo), (17) 

U+ = U- = f 2 (U ,N ), (18) 

M+=M- = f 3 (U ,Mo), (19) 

M- =M+ = U{Uq,Mq), (20) 



with the functions on the right-hand sides depending on 
Uq, Mq and other model parameters, the corresponding 
QSH and LAF solutions can be immediately determined 
as well. In particular, the QSH solution is given by 



A+ = -A_ = fi(U ,Mo), (21) 

U+ = U. = f a (U ,Afo), (22) 

M + =M- = f 3 (U ,Mo), (23) 

M+ = -M- = U(U ,M ), (24) 

whereas, the LAF solution is given by 

A+ = A_ = 0, (25) 



u± = /2(^o,M))T/i(t/o,M)) | (2g) 

MU ,Mo)±U(UoM , 07 , 
A± = ^ ' ( 27 > 

M-=M+ = 0. (28) 

IV. RESULTS 

In this section, we present our main results. We start 
from the analysis of the gap equations in the case of van- 
ishing top-bottom gates voltage imbalance, Uq = 0. In 
particular, we reveal the emergence of the hybrid state at 
nonzero bare nematic parameter Mq and describe how it 
gradually turns into a pure nematic state with increasing 
Mq. We then discuss the complete phase diagram in the 
Uq-Mq plane. 

A. Nematic and spin-valley symmetry-broken 
states at Uo — 

In the hybrid QSH, QAH, and LAF states (with an ad- 
dition of nematic ordering), we find that in general the 
nematic parameters depend on spin and/or valley indices, 
i.e., A+ 7^ A"- and M s ^ 0. However, in the case of the 
vanishing top-bottom gates voltage imbalance, Uq = 0, 
the structure of w^ s simplifies: M+ = M- and M s = 0. 
The latter follows from the structure of Eqs. (|12[) . (|T3"|) . 
and (fTS"]) . In order to clarify the role of a bare nematic 
parameter Ao in competition between the nematic and 
the gapped states, we first study this case. Then, the 
numerical analysis of gap equations (|T2"j) and ([T3| shows 
that solutions of the two types are possible. One of them 
is the nematic solution with only w^ s 0, and the other 
is a hybrid solution with both A^ s and w^ s nonzero. Note 
that a gapped solution with only A^ s ^ is impossible 
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FIG. 1: (Color online) The order parameters A s and M a as 
functions of the bare nematic term A/o for the nematic solu- 
tion (solid line) and for the hybrid QSH solution (dashed and 
dotted lines) at Uq = 0. 

when a nonzero bare nematic parameter Ao is present, 
therefore, gap and nematic parameters coexist in the hy- 
brid solution. 

At vanishing Uq and Ao , the Q VH solution has a higher 
energy than the QAH, LAF, and QSH solutions. The ad- 
ditional energy cost is due to the Hartree term in 
the gap equation (fT2j). This term is negative because it 
reflects the energy price associated with the electric field 
between the layers of graphene in the QVH phase. Thus, 
there arc two main solutions of the gap equations (fT2")) 
and (|13[) . the pure nematic and one of the hybrid solu- 
tions, QAH, LAF, or QSH. The latter are degenerate in 
energy. Instead of discussing all the degenerate solutions, 
in the rest of the paper, we will concentrate only on the 
hybrid QSH one. 

In the pure nematic solution, the nematic order param- 
eter A" s does not depend on the valley or the spin. Its 
dependence on the bare parameter Ao is shown in Fig. [T] 
At small Ao, it can be approximated by a linear depen- 
dence A± = A'ofi + bMo with the offset M s =4.61 meV 
and the slope b = 3.18. For the hybrid QSH solution, the 
gaps A s and M s are also plotted in Fig. [T] We see that 
the bare nematic parameter inhibits the gap. 

Comparing the energy densities of the nematic and 
hybrid QSH solutions, we find that the hybrid solution 
always has a lower free energy whenever it exists, see 
Fig. [21 From Fig. [TJ we see that the dynamical nematic 
order parameter M s for the hybrid QSH solution van- 
ishes as Ao — > 0. Therefore, the ground state of unbi- 
ased bilayer graphene at Ao = Uq = is a spin- valley 
symmetry-broken gapped state. Whereas, the dynami- 
cal nematic order parameter M s grows, the gap A s grad- 
ually decreases with increasing the bare nematic param- 
eter Ao- Eventually, the hybrid QSH solution smoothly 
turns into the pure nematic solution at the critical point 
Mq t ps 1.26 meV where the gap A s turns to zero. Since 
the QSH state breaks the spin SU (2) symmetry, whereas, 
the nematic state does not, we conclude that the cor- 
responding transition is a second-order phase transition 
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FIG. 2: (Color online) The energy density as a function of the 
bare nematic term Ao for the nematic solution (solid line) and 
the hybrid QSH solution (dashed line) at Uo — 0. 
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FIG. 3: The phase diagram in the Ao-Uo plane. 



(and not a smooth crossover). 



B. Phase diagram 



FIG. 4: (Color online) The order parameters as functions 
of the bare nematic term Ao for the hybrid QSH solution 
(thick red lines) and hybrid QVH solution (thin blue lines) at 
Uo = 0.4 meV. 

In the subcritical region (i.e., at small values of Ao 
and Uo) of the phase diagram at a fixed Uo 7^ 0, in- 
creasing the bare nematic term A/"o induces a substantial 
dynamical nematic order parameter w^ s , see Fig. |4] for 
Uo = 0.4 meV. Consequently, in this case, the hybrid 
QSH state has spin-symmetric order parameters U s and 
A/" s , and spin asymmetric order parameters A s and A4 S . 
In the supercritical region, the hybrid QVH state with 
spin singlet U s and Af s parameters (see Fig. @| has the 
lowest energy. As expected, the hybrid QVH state be- 
comes almost a pure nematic state at small Uo and large 
A/o and almost a pure QVH state at large Uo and small 

The critical line in Fig. [3] separates the hybrid QSH 
and QVH phases. Its top part corresponds to a first- 
order phase transition (the gaps of the hybrid QSH and 
QVH states are discontinuous across the critical line), 
which ends in a critical end point at A/q* ~ 0.36 meV and 
Uq ps 0.95 meV. The rest of the critical line corresponds 
to a second-order phase transition. 



Using a coupled set of gap equations for the nematic 
and several types of spin-valley symmetry-breaking or- 
der parameters in bilayer graphene, we can now obtain a 
phase diagram in the plane of two parameters: a strain- 
induced bare nematic term Ao and a top-bottom gates 
voltage imbalance Uo- This is performed by sweeping 
through the corresponding two-dimensional space of the 
bare parameters, comparing the energies of all solutions, 
and determining the ground state at each point. The 
results are summarized in Fig. [3J 

To understand the overall topology of the phase dia- 
gram, it is useful to first look at the competition of the 
non-nematic phases at A/o = 0*£~— In this case, increas- 
ing the voltage imbalance Uo has a tendency to suppress 
the QSH state gap A s and to induce an increasing ad- 
mixture of the U s gap in this state. At a sufficiently large 
Uq, a pure QVH state with a spin singlet U s gap takes 
over as the ground state. 



V. SUMMARY 

The current paper shows that bilayer graphene has a 
rich phase diagram in the plane of the bare nematic term 
Ao and the voltage imbalance between the layers Uo~- 
It is argued that the ground state should be one of sev- 
eral hybrid states with spontaneous spin- valley symmetry 
breaking. 

We find, in particular, that one of the hybrid states, 
QAH, LAF, or QSH state is the ground state in the sub- 
critical region (i.e., small values of Ao and Uo), and that a 
hybrid QVH state is the ground state in the supercritical 
region (i.e., large values of Ao and/or Uo). The transi- 
tion between the two phases is either a first-order phase 
transition (at small Ao and large Uo), or a second-order 
phase transition. The two parts of the critical line meet 



at a critical end point, which is an interesting prediction 
in its own right. 

Note that, in the model with the Coulomb interaction 
used in this paper, three of the proposed hybrid states 
(i.e., QAH, LAF, and QSH) appear to be degenerate in 
energy. This degeneracy is expected to be lifted after 
the inclusion of additional symmetry-breaking contact in- 
teractions. Indeed, as follows from the renormalization- 
group studies in the normal phase^^ ' 22 ' 23 there are sev- 
eral types of such interactions that may become relevant 
for the low-energy dynamics. The corresponding general- 
ization of the phase diagram is of great interest and will 
be addressed in the future^ 

Based on the general arguments reviewed in Ref. [2^, 
one should expect that, although gapless edge states are 
present in (hybrid) QAH and QSH states, there arc no 
such states in (hybrid) QVH and LAF ones. The role 
of the edge states in present dynamics will be considered 
elsewhere. 

It is interesting to note that certain two-dimensional 
fcrmionic models with a quadratic band-crossing point 
similar to bilayer graphene, but with a <i-wave symmetry, 
can also have a hybrid phase with coexisting spin- valley 
symmetry-breaking and nematic order parameters^ 

The experimental data in Ref. [Ill were fitted with the 
nematic parameter on the order of 6 meV, and it was 
argued that the required strain to explain the observed 
effects should be too large. According to our results, in- 
teractions strongly enhance the bare nematic parameter 
(e.g., ~ 6.3 for Af Q ~ A^ r at U a = 0). There- 

fore, even a rather small strain can lead to the nematic 
state. This suggests that different bilayer samples in ex- 
periments may have very different physical properties, 



depending on the specific conditions under which they 
were created. 
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Appendix A: Energy density 

According to Eq. (2.9) in Ref. [13, the energy den- 
sity of a three-dimensional electron gas can be calculated 
through an integral of the electron Green's function. It 
is not difficult to check that, with obvious modifications, 
this formula is valid also for electron quasiparticles in bi- 
layer graphene. Alternatively, the energy density can be 
obtained by evaluating the Baym-Kadanoff functional on 
the solutions of the gap equations. Then, we have 



die 
2^ 



d 2 k 



tr [ (u + £[/ot 3 + AVi + H ) G ?s (w, k) ] 



(Al) 



where Hq is the free Hamiltonian of bilayer graphene. Using the expression for the propagator in Eq. (|10p . we derive 
an explicit expression for the energy density, 



dii) 
2^ 



d 2 k c 2 + U A £s + ^ + k2(Afo+w J ^ )cos2 ^ + M 



„2 _ J±_ 

U £s Am 2 



k 2 W£ s COS 1^>h 



iO 



(A2) 



which can be used for any state with dynamically generated gaps and nematic parameters. Note that the integral 
over u diverges in this expression for £. Therefore, we should subtract the "vacuum" energy of the trivial solution 
with A^s = w^ s = 0. Then, we obtain 



Ef du> 
/ 9^ 



d 2 k 



2tt J (2tt) 



UoA. 



fc 2 (7Vo+m; a ) cos 2ip k 
2m 



fc 4 

4 m 2 



; 2 - a 2 s 



v 2 - -K, 



k 2 W£ s COS 2t^fc 



iO 



1 k 4 

Am* 



id 



(A3) 



Performing the Wick rotation and integrating over u>, we arrive at the expression, 
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The integral over the momentum in this expression is is determined by the range of validity of the effective 
logarithmically divergent at large k. In the calculations, Hamiltonian (TTJ) (recall that A = 7i/4). 
therefore, we use a finite cutoff at fc max = \/2mA, which 
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